## Authors: Alexander Herzog and Kenneth Benoit
## Date: May 30, 2015
## Replication file for JOP article "The Most Unkindest Cuts: Speaker Selection and Expressed Government Dissent During Economic Crisis"

rm(list = ls(all = TRUE))

library(rjags)
set.seed(42)

### PATHS ########################################
dataSub <- "./generated_data/working_data_speakers.RData"
##################################################

# MCMC values
mcmc.adapt <- 500
mcmc.burn <- 1000
mcmc.sample <- 2000
mcmc.thin <- 2
mcmc.chains <- 2

# JAGS model 
JagsModel <- "model{
	for (i in 1:N)
	{
		y[i] ~ dnorm(mu[i], tau)
		mu[i] <- a[m[i]] + inprod(b, X[i,1:K])
	}
	
	for (j in 1:J)
	{
		a[j] ~ dnorm(mu.a, tau.a)
	}

	for (k in 1:K)
	{
		b[k] ~ dnorm(0, 0.0001)
	}

        tau <- pow(sigma, -2)
        sigma ~ dunif(0, 1000)

        mu.a ~ dnorm(0, 0.0001)
	tau.a <- pow(sigma.a, -2)
	sigma.a ~ dunif(0, 1000)


}"


# function for JAGS data for multilevel models
forJagsData <- function(DV, xMatrix, groupID) {
    list(
        y = DV,
        X = xMatrix,
        N = nrow(xMatrix), 
        K = ncol(xMatrix),
        J = length(unique(groupID)),
        m = groupID
    )
}

load(dataSub)
data <- dataSub


# Model 1 without interaction effects
# -----------------------------------
print("Estimating outcome model 1")

eq <- textscore ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis
X <- model.matrix(eq, data)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(data$textscore, X, data$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos1.posterior <- coda.samples(model,
                                 variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(pos1.posterior,file="./estimated_models/outcome_model1_posterior.RData")


# Model 2 with interaction effects
# ---------------------------------
print("Estimating outcome model 2")

eq <- textscore ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*unemployment + crisis*safetyScaled

X <- model.matrix(eq, data)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(data$textscore, X, data$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos2.posterior <- coda.samples(model,
                                 variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(pos2.posterior,file="./estimated_models/outcome_model2_posterior.RData")


# Model 3 with three-way interaction effects
# ------------------------------------------
print("Estimating outcome model 3")

eq <- textscore ~ 0 + unemployment + safetyScaled + seniorityYearsScaled + government + cabMember + crisis + crisis*unemployment + crisis*safetyScaled + crisis*government + government*unemployment + government*safetyScaled + crisis*government*unemployment + crisis*government*safetyScaled

X <- model.matrix(eq, data)

model <- jags.model(
    textConnection(JagsModel),
    data = forJagsData(data$textscore, X, data$m),
    n.chains = mcmc.chains,
    n.adapt = mcmc.adapt)

update(model, mcmc.burn)

pos3.posterior <- coda.samples(model,
                                 variable.names=c("a", "b", "sigma", "mu.a", "sigma.a"),
                                 n.iter = mcmc.sample,
                                 thin = mcmc.thin)

save(pos3.posterior,file="./estimated_models/outcome_model3_posterior.RData")

